!! Copyright (C) 2009,2010,2011,2012  Carlo de Falco
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Carlo de Falco                    <carlo.defalco@polimi.it>

program main   
  use mod_petscintf, only: &
       mod_petscintf_constructor, &
       petsc_init, &
       petsc_solve, &
       petsc_print, &
       petsc_cleanup, &
       mod_petscintf_destructor

  implicit none

  
  integer :: &
       ir(2998), ic(2998), &
       gn(1000), ii, jj

  real :: &
       a(2998), sol(1000), &
       rhs(1000)


  !! init
  call mod_petscintf_constructor ()
  call petsc_init (1000)

  !! create matrix, solution and rhs
  do ii=1,1000
     ir(ii) = ii-1
     ic(ii) = ii-1
     gn(ii) = ii-1
     a(ii)  = 2.0
     
     sol (ii) = 1.0
     rhs (ii) = 1.0
  end do

  do ii=1,999

     ir(ii+1000) = ii-1
     ic(ii+1000) = ii
     a(ii+1000)  = -1.0

     ir(ii+1999) = ii-1
     ic(ii+1999) = ii
     a(ii+1999)  = -1.0

  end do

  call petsc_solve (1000, ir, ic, a, gn, sol, rhs)

  !! print matrix, rhs and solution for debugging
  call petsc_print ()

  !! free matrix, solution and rhs
  call petsc_cleanup ()

  !! quit
  call mod_petscintf_destructor ()

end program main

